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We study time-optimal protocols for controlling quantum systems which show several avoided 
level crossings in their energy spectrum. The structure of the spectrum allows us to generate a 
robust guess which is time-optimal at each crossing. We correct the held applying optimal control 
techniques in order to find the minimal evolution or quantum speed limit (QSL) time. We investigate 
its dependence as a function of the system parameters and show that it gets proportionally smaller 
to the well-known two-level case as the dimension of the system increases. Working at the QSL, 
we study the control fields derived from the optimization procedure, and show that they present a 
very simple shape, which can be described by a few parameters. Based on this result, we propose a 
simple expression for the control held, and show that the full time-evolution of the control problem 
can be analytically solved. 


I. INTRODUCTION 

The development of future communication and infor¬ 
mation technologies is expected to rely heavily on the 
precise manipulation of physical systems at the nano- 
and sub-nanoscale. For this reason, coherent control of 
quantum systems has become a major goal in physical 
sciences over the past decades. In this context, the 
design and implementation of quantum control methods 
has raised a lot of interest, and many theoretical m 
and experimental HE] works have been devoted to this 
subject. 

Controlled quantum operations are tipically intended 
to be performed in the fastest possible way, in order to 
avoid unwanted environmental effects which can destroy 
the coherence properties of the system. Because of 
this, deriving time-optimal control protocols is a goal of 
major importance. This task is usually tackled by means 
of quantum optimal control (QOC) theory. There, the 
typical problem is to derive the shape of the control 
field A (f) required in order to optimize a particular 
dynamical process for a quantum system described by a 
Hamiltonian H( A). For example, a typical objective in 
quantum control is to perform a transition from a given 
initial state \ipo) to another goal state \i/j g ). In some 
cases, the optimization can be carried out analitically 
EH, but most generally its approached numerically 

mini. 

One of the weak points of the usual algorithms 
employed for solving QOC problems (such as Krotov or 
GRAPE [131 3, although interesting alternatives have 
been proposed recently is, m) is that the solution 
for the field A (t) often appears to be hardly realizable 
in practice m • This originates from the fact that the 
value of the held at each instant acts as an independent 
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control (i.e. there are no constraints derived from the 
truncation of a given basis set of functions iTS]). Of 
course, this allows for a much faster convergence of 
the optimization procedure, but the resulting control 
held can present non-smooth fluctuations which would 
require a large held band with to be implemented. 
Moreover, from a theoretical perspective, the complex 
shape of the held usually prevents us to understand the 
physical mechanisms involved in the control processes. 
Nevermind this feature, QOC theory has been proven 
to show deep connections with the fundamentals of 
Quantum Mechanics. Caneva et al. m studied the 
performance of QOC in various systems as a function of 
the (hxed) evolution time T that is fed to the algorithm. 
They found that the optimization converged succesfully 
only when T was above certain threshold, which they 
identihed with the quantum speed limit (QSL) time, 
Tqsl ■ The concept of QSL was introduced originally 
by Mandelstam and Tamm [20]) who showed that a 
generalization of the usual time-energy uncertainty 
relation imposed bounds on the speed of evolution of 
a quantum system. Since then, many authors have 
explored the QSL in various situations [2TH32] . 

In this work we study time-optimal control processes 
in quantum systems which show several local two-level 
interactions in the form of avoided crossings (ACs) in 
their energy spectrum. This situation is of interest in 
a wide variety of quantum mechanical systems, such 
as molecular dipoles interacting with electric fields 
[331133] , ultracold atoms in optical lattices [35], Rydberg 
atoms [36] and superconducting qubits [37]. Using this 
particular interaction between the states of the system, 
we generate initial guesses for the control protocols 
using piecewise-constant functions derived from previous 
studies 138;, [39]. In a recent work, we studied the QSL 
time for these protocols by using QOC, and showed 
that the calculated QSL time is in general smaller than 
the sum of the optimal times for each avoided crossing 
[40] . Here, we focus on the analysis of the optimal 
control protocols which lead to such speed-up. For that 
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purpose, we investigate numerically the control fields 
that generate the time-optimal evolution and find that 
they can be fully characterized by just a few parameters. 
This allows us to propose a simple analytical dependence 
for the control field. Finally, for this protocol we show 
that the full time-evolution can be analytically solved, 
and the results are in full agreement with the numerical 
optimization. 

This article is organized as follows. In Sec. [Il] we 
present the model of a quantum system showing an 
avoided crossing (AC) in its energy spectrum and de¬ 
scribe its most importante features. We then expand 
this minimal model to include several ACs, for which 
we present the actual model Hamiltonian of our inter¬ 
est. We also discuss the control problems that can be 
posed for this system, and present an intuitive solution. 
In Sec. |III| we present the basics of optimal control the¬ 
ory and describe its implementation in quantum systems. 
In Sec. |IV| we briefly discuss the results obtained by im¬ 
plementing QOC for control processes involving several 
ACs, and the study the QSL time as a function of the 
different parameters of the system. In Sec. [V] we de¬ 
scribe in detail the optimal control fields we obtain, and 
discuss the physical mechanisms involved in the observed 
speedup. Based on this analysis, we propose a simple an¬ 
alytical expression for the optimal control field, and show 
that the associated Schrddinger equation can be solved 
analytically. Finally, Sec. |VI| contains some concluding 
remarks. 


II. MODEL, AVOIDED CROSSINGS AND 
CONTROL PROTOCOLS 

In this section we present the models which describe 
the systems of our interest, which show avoided level 
crossings in their energy spectrum. We propose simple 
control protocols for achieving state transfer and discuss 
its time-optimality. 

A. A single avoided crossing 

We first consider a quantum two-level system described 
by the following Hamiltonian matrix 

H W = + Act 2 = ^ a q)> (!) 

which is written in the basis {10) , 11)}. These states are 
usually called the diabatic states of the system, which 
diagonalize the Hamiltonian when the control parameter 
A -A ±oo. In general, the eigenvalues {Ek} (k = 0,1) of 
H form a hyperbolae in the (A, E) plane, whose vertex 
represents an avoided crossing (AC) with an energy gap 
A. This spectrum is depicted in Fig. [l] (a). The eigen¬ 
states of H as a function of A form the adiabatic basis 


{|<?a) , |ca) } and have an asymptotic correspondence 
with their diabatic counterparts, i.e., \g~oo) = |0) and 
|e_oo) = |1) (and viceversa for A -A ± 00 ). 

The model presented above is ubiquitous in quan¬ 
tum mechanics as it accounts for many interesting 
phenomena, such as Landau-Zener transitions 02: 
Landau- Zener- Stuckelberg interferometry [42) and quan¬ 
tum phase transitions II. 3 1 . We are interested in the con¬ 
trol problems that can be formulated when A is regarded 
as a fixed parameter, and A can vary in time. A famous 
example is the problem of driving this system from |g+A 0 ) 
to Iff-Ao) i n the shortest possible time, for some A 0 € R. 
Interesting discussions about the solution to this prob¬ 
lem, which include numerical, experimental and analyti¬ 
cal studies can be found in the literature [21 [3H [44H46]. 
Here we will focus in a particular result. When Ao —> 00 , 
the control problem stated above reduces to the full pop¬ 
ulation transfer between |0) between |1). The optimal 
time for such process is given by 


and can be achieved simply by setting A (t) = 0 from 
t = 0 to t = Tg X \ given of course that |^>o) = |0). Then, 
the state can be frozen in the final state for t > by 
applying a quench from A = 0 to some value |A| A. 
An example of this type of control field is depicted in Fig. 

0(b)- 

B. Multiple avoided crossings 

The two-level model described in the previous section 
can be extended and generalized to account for the pres¬ 
ence of several ACs in a many-level scenario. Here we 
construct a model for such situation. Consider an N- 
level system with the following Hamiltonian 

[^] 

H n { A) = ^2 (A-ne 0 )|2n)(2n| 

n— 0 

+ ^ ^ Tl € 0 1 2 t 7 , ~\~ 1 ) ( 2 ? 7 - + 11 

n—0 
N ~ 2 A 

± ^ (|n)(n± t| ± |n± tXnl), (3) 

n— 0 

where [x] denotes the integer part of x and {|n}} 
is the basis of diabatic states. When A n = 0 for 
n = 0,1,..., N — 2, the Hamiltonian is diagonal in that 
basis, and the energy spectrum consists merely on a 
series of horizontal and diagonal straight branches with 
degeneracies at values of A ij = (i + j)e 0 corresponding 
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Control parameter A 


Figure 1. (a) Energy spectrum for the two-level Hamilto¬ 

nian (JTJ as a function of parameter A. (b) Control field as a 
function of time for a simple realization of complete popula¬ 
tion transfer between the diabatic states, (c) and (d) same 
as (a) and (b) but for the three-level Hamiltonian 0 - (e) 
Schematical representation of the energy spectrum of Hamil¬ 
tonian Hn{ A), c.f. Eq. 0- as a function of control parameter 
A. In the most general setting, the spectrum shows N — 1 
avoided crossings separated by eo, each of which generate a 
coupling between states | n) and |n + 1) with n < N — 1. 


to states 1 2i) and |2 j + 1). If one of the couplings is non¬ 
zero, say A ra 7^ 0, the degeneracy at X n = ne o is lifted 
and an AC is generated with a minimum energy gap of 
A„. As a consequence, transitions between the states | n) 
and | n + 1) become permitted. The overall shape of the 
energy spectrum for this model is schematically depicted 
in Fig. |T] (e). Note that, when all the interaction rates 
are non-zero, the number of ACs equals N — 1 . In a 
regime where eo A„ for all n , this model is very 
convenient for analyzing dynamical processes which are 
dictated by local two-level interactions. This can be seen 
as follows: if the system is initially preparred in some 
state | n) and the control parameter A does not deviate 
much from the position of the corresponding AC (i.e. 

| A — A„| -C A n [50]), then the dynamics of the system 
is effectively confined to a two-dimensional subspace, 
as the remaining N — 2 levels can be adiabatically 
eliminated m ■ This is the key characteristic of our 
model, and we will expand on its consequences later on. 

Evaluating Eq. 0 for N = 2 we recover the two- 
level (one AC) Hamiltonian of Eq. 0 Taking the next 
step in complexity, the case N = 3 renders the following 


Hamiltonian matrix 

/A^ 0 \ 

H 3 (X) = t 0 t > ( 4 ) 

V 0 t X~e 0 J 

which has two ACs, one at Ao = 0 and other at 
Ai = eo- The corresponding gaps are Ao and Ai 
when eo Ao,Ai. The energy spectrum for this 
case is depicted in Fig. [I] (c). This model has been 
widely studied in many different contexts pTTTFdTJj . as 
it is suitable for describing a three-level atom in a A 
configuration. Note that, in that case, the parameters 
Ao and Ai correspond to detunings between the energy 
levels and the frequencies of two external laser fields, 
which are generally regarded as the control fields, while 
A and eo are related to the bare energy splittings. In this 
work this is not the case, as the off-diagonal couplings 
are fixed and we implement control protocols by variying 
solely A (t). 

For this multiple AC model, we are interested in con¬ 
trol processes which connect diabatic states of the sys¬ 
tem. Without loss of generality, we consider the ini¬ 
tial state \ipo) = |0) and define the process Pk as the 
one which drives the system to the state | K), with 
0 < K < N — 1 (generalization to a different dia¬ 
batic initial state is straightforward). Our goal will be 
to find the control held A x{t) which generates Pk in 
a time T. Note that, if the ACs are sufficiently iso¬ 
lated, a solution exists which is independent of N. Based 
on the discussion above, an straightforward, yet power¬ 
ful solution can be drawn [39]. By succesively setting 
A (f) = X n constant during time intervals of length tt/ A n , 
with 0 < n < K — 1 the dynamics navigates through 
the I\ ACs turning them on and ensuring full popula¬ 
tion transfer one at the time. The system then evolves 
through the sequence |0) —► |1) —>■... —► | K). Note that 
the shape of the control function is then characterized 
by a series of sudden changes of the value of A, an so 
we name it a “sudden switch” field, X^\t), which is de¬ 
picted for K = 2 in Fig. [l] (d). The total evolution time 
for this protocol equals 



K -1 


E 

n =0 


7r 



(5) 


A number of observations are in place. First, note that 
the protocol proposed is not unique, since the process 
can also be realized by adiabatically changing A (t) as 
to navigate through the ACs. The system then also 
evolves sequentally between diabatic states, but much 
more slowly [3SJ. Also, is important to point out that 
we have constructed the model in Eq. 0 in such a way 
that the degeneracies between states | n) and | n + k) (for 
k 1) are exact, and cannot be lifted. For this protocol, 
this means that there is only one path in the energy 
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spectrum between |0) and | K), which involves exactly 
K ACs. We point out that we do not lose generality by 
making this assumption: if there were a shorter path 
between those states, it would be equivalent to a process 
Pl with L < K, which is accounted for in our model. 
Finally, let us remark that the total protocol time given 
by Eq. © is a sum which terms are of the form of Eq. 
([2]), an so we can state that the control saturates the 
QSL bound at each AC. In the following, we investigate 
wether this feature implies that the whole protocol is 
itself optimal or not. 


III. OPTIMAL CONTROL THEORY 


cost of the control process mm- Note that the factor 
a(t) allows for selective weighting at different times, thus 
allowing the induction of certain special features in the 
control field (e.g. a given shape, or its initial and final 
values EZI) As a consequence, we also wish to maximize 

J 2 [A] = — f a(t)\ 2 {t)dt. (9) 

Jo 

Finally, a restriction has to be imposed to the joint 
maximization of Eqs. © and (pl), in order to guarantee 
that the dynamical equation (|6J) is satisfied at all times. 
For that purpose we introduce an auxiliary state |y(t)) 
as a Lagrange multiplier so that we seek to maximize a 
third functional 


Here we outline the theoretical formulation of a basic 
QOC problem, its solution and a feasible method for its 
numerical implementation. Details on this derivation 
can be found for example in Ref. m- 

Consider a quantum system described by a state | if>(t)) 
in a finite-dimensional Hilbert space H undergoing an 
evolution dictated by the Schrodinger equation (we take 
h = 1 from here on) 


J 3 bcVh A] = -21m | J 


(xWI 




In conclusion, bringing together expressions ([8]) 
through (10) we get that the QOC problem is casted 
as the maximization of the functional 


I^W) = #(*) I^W) > ( 6 ) 

satisfying |i/>(0)) = |^o)- We supose that the Hamiltonian 
of the system H (t) has the form 

H(t) =H 0 + A (7) 

where Hq and Hi are the free (or drift) and interaction (or 
control) Hamiltonians, which are time-independent, and 
we define real-valued function A (t) as the control field. 
The general QOC problem is formulated as follows: given 
H 0 , Hi, an initial state |^>o)> an hermitic operator P and 
a total evolution time T, we wish to find A(f) such that 
the system, initially prepared in \i/jq), evolves to a state 
| tp(T)), in which the expectation value of P is maximal. 
Formally, we wish to maximize the following functional: 


A] = Ji[V>] + Ja[A] + J 3 [x,i>i A]. (11) 

Optimization of this functional is achieved by imposing 
SJ = 0, which renders three independent equations (one 
for each variable of the functional). First, solving S y J = 
0 trivially gives Eq. as expected from the inclusion 
of the Lagrange multiplier. Then, the relation 6^>J = 0 
takes us to the following equation for the auxiliary state 

lx) 

4 IX(*)> = H (*) I X(t)) and |x(T)> = P \^(T)) . (12) 

Note that this expression is the Schrodinger equation 
for state |x(i)), with boundary condition given by its 
final value, |x(T)}. Finally, by solving 5J\ = 0 we an 
expression can be derived for the control field 


J^\ = ^{T)\P\^{T)) (8) 


A(*) = Z 77 \ lm {(x(t)\H i \'ip(t))} . (13) 

a(t) 


In the following we will restrict ourselves to the case 
in which the goal of the QOC problem is to maximize 
the probability of reaching a certain goal state | ip g }- For 
that purpose, the operator P is defined as the projector 
P = Hg){tp g \ and so Ji [ijj\ = \(ip(T)\i{; g )\ 2 . 

In order to correctly formulate the QOC problem, two 
additional conditions have to be imposed. The first one 
is the minimization of the quantity f 0 a(t)A 2 (t)dt, where 
a(t ) is a weight function. This requirement is essential 
in order to prevent the divergence of the total energy 


The problem of obtaining a set {\ip(t)), |x(^)) > A(t)} 
that simultaneously solve Eq. © together with Eqs. 
(12) and (13) is, of course, impossible to tackle an¬ 
alytically. Instead, an iterative algorithm has to be 
implemented. Here we briefly describe a widely used 
method [52J [S3], due originally to Krotov pT3]: (i) the 
procedure starts by choosing an initial guess Ao (t) 
for the control field; (ii) using that field, the initial 
state |^o) is evolved according to Eq. © from t = 0 
to f = T; (iii) the boundary condition for |x(t)) is 
set by projecting |x(T)) = P\ip(T)), and the state 
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is evolved backwards also following Eq. |6]), from 
t = T to t = 0; (iv) the state \i/jq) is now propagated 
forward again, but the field is updated following the rule 
A (t) —> A (t) + ^ylm{(x(t)| Hi IV’(i)}} at each instant; 
(v) steps (iii) and (iv) are repeated N times until a 
certain threshold is reached for the value of the cost 
functional J\. 


IV. QSL IN A SYSTEM WITH MULTIPLE ACS 

In this section we numerically investigate the QSL time 
for the control processes described in the previous sec¬ 
tion. For that purpose we use optimal control techniques, 
inspired by the basic idea introduced by Caneva et al. 
m that the optimization performance is limited by the 
maximum speed allowed by quantum evolution. The ba¬ 
sic procedure is as follows. First, we fix the state dimen¬ 
sion N and choose a control process Pk for the model 
described in Section |TT| Then, we run the optimization 
algorithm in order to find the control field Xx(t) which 
generates the desired process, for different values of the 
total evolution time T. In each run, this procedure takes 
as an input the value of T and an initial guess for the 
field A® (t). In order to choose these inputs, we take ad¬ 
vantage of the physical features of the model discussed 
in the previous section. The values of T were taken from 
an interval centered around Tg K \ cf. Eq. Note 

that, if the ACs are well isolated, we are certain that the 
sudden switch field generates the desired process when 
T = Tg K 1 . Similarly, the initial guess for the control 
function were chosen to be close to the sudden switch 
field. Actually, we used 

(t.) = a (t) A^ S) (b t) + c(t) (14) 

where b is a parameter which shrinks or expands the 
shape of the function to fit the total evolution time (i.e. 
b = 1 when T = Tg K •*), while a{t) is a function which 

smooths the discontinuities of X y K and c(t) is a small 
linear correction. The latter functions are introduced 
in order to force the algorithm to take a minimum 
number of steps (of the order of 100) before the required 
convergence is achieved. 

Each run of the algorithm finishes after a fixed number 
of steps, or when the process is sufficiently converged. As 
discussed in Section [Till this is determined by evaluating 
the value of the infidelity at each step to, which is defined 
as 


T = 1 — 


v> s iv> (m) cn )| 2 = i-(is) 


where J\ is the functional of Eq. 0, and 

is the state of the system obtained at step m of the 



Figure 2. (color online) A typical procedure for determining 
the QSL time Tqsi for a particular control process. This 
example corresponds to a Pk= 2 process in a Hilbert space of 
dimension N = 3, cf. Eq. 0 , where eo = 10Ao and Ai = Ao. 
(a) Infidelity T, n as a function of the step number m for each 
optimization, (b) Second derivative 1'^ of the curves in (a). 
The dotted line corresponds to the minimum value of T which 
asymptotically renders T'Xn < 0 and so is identified as Tqsl 
within error margin. Thin full lines correspond to T < Tqsl 
while, thick full lines to T > Tqsl- 


algorithm. The function I m decreases monotonically as 
to increases, but its shape and asymptotic behaviour 
depends critically on the input parameters. In Fig. [2] 

(a) we plot this function for a particular case, as an 
example. We argue, as in Ref. [IS] that the infidelity 
cannot decrease indefinitely if T is smaller than the QSL 
time. In that case, I m should look asymptotically flat. 
We use this feature to obtain the estimator of the QSL 
time Tqs L . Formally, for each value of T we look at the 
second derivative of X m (with respect to to), see Fig. [2] 

(b) , and analyze its sign. Then, the minimum value of 
T which gives X"{k) < 0 asymptotically, is chosen as the 
QSL time. 

We now turn our focus to the model of Eq. © 
which presents two ACs. As discussed in the previous 
Section, this is the next step in complexity following 
the analytically solvable two-level system. We begin 
by considering the QSL time for process Pi, for which 
the system starts in state |0) and evolves to 11), in the 
minimum possible time. Note that this process involves 
just one AC, as seen from the sudden-switch protocol 
introduced in the previous section. In Fig. [3j we plot the 
calculated QSL time Tq^ l for this case as a function of 
6q , the parameter which measures the distance between 
the ACs in the energy spectrum (see Fig. [I]), for fixed 
values of interaction parameters Ao,Ai. There, it can 
be seen that Tqq L is larger than Tg 1 ^ = 7r/Ao for small 
values of eo- This is reasonable in this regime, since 
the ACs interact considerably, which leads to significant 
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Figure 3. (color online) (a) QSL time calculated from the 
optimal control procedure (see text for details) as a function 
of eo/ A o, for processes Pi (crossing one AC) and P2 (crossing 
two ACs). Dashed lines correspond to expression ©. i-e- 
the time required by the sudden-switch protocol in each case, 
Tg 1 ^ and (b) Ratio between the calculated QSL time 

Tqsl for processes Pk and the corresponding sudden-switch 
protocol evolution time Tg K ■* as a function of K. The dashed 
line show the A' -1 scaling of the data points. For all cases, the 
distance between the ACs was set to eo = 10A, and A„ = 1 
for all n = 0,1,..., K — 1. 

variations of the interaction rates (see Ref. [47] for more 
details). Away from that regime, Tqs l converges to 

Tg , which is the well-known result for the two-level 
system. This is a sound result, since only the states 
|0) and |1) are involved in the process. However, it is 
interesting to point out that this behaviour allows us to 
quantitatively define the regime in which the ACs are 
well isolated. In the case shown in the figure, for which 
A 0 = Ai, this is achieved for eo/A a 5. 

Next, we discuss control process P2, which involves 
both ACs. Following the same procedure as for the 
previous case, we get the results of Fig. [ 3 ] (a). There, 

( 2 ) 

it can be seen that the estimated QSL time Tgg L is 

smaller than the sudden switch evolution time Tg . 
Remarkably, this result holds in all cases, even for large 
eo- The difference between Tgg L and our prediction 
is larger for small eo, and decreases as the ACs are 
brought apart. However, for e 0 /A 0 as large as 100, 
the difference is still larger than 7%. This striking 


behaviour indicates that the QOC optimization can 
generate successful (i.e., with arbitrary fidelity) control 
processes which are significantly shorter in time than the 
double sudden-switch, a process wich is time-optimal at 
each AC, as discussed above. We point out that this 
behaviour persists even when the relative magnitude 
of the gap sizes is modified [J0|. We will analyze the 
physical mechanisms that cause this speed-up in the 
next section. 

Finally, we address the results obtained for the QSL 
time for control processes involving more than two ACs, 
i.e. Pk with K > 2. Applying the same procedure out¬ 
lined in the previous paragraphs, we obtained Tqsl for 
various values of the number of avoided crossings K in¬ 
volved in the process. In Fig. [3] (b) we plot the ratio 

between Tqsl and Tg K ' > as a function of K. There, it 
can be seen that the optimal evolution time (measured 
with respect to the corresponding sudden switch protocol 
evolution time) decreases as the number of ACs involved 
increases. This means that, as more ACs get involved in 
the evolution, the connection between diabatic states can 
be performed faster. However, the improvement reaches 
a saturation point for large values of K. 

V. OPTIMAL CONTROL FIELDS 
A. Numerical analysis of the optimized fields 

We now turn to analyze the shape of the control fields 
derived via the optimization procedure outlined in the 
previous section. We will focus on the optimal fields 
obtained for T = Tqsl , but for larger evolution times its 
description is similar. In Fig. [4] (a) and (b) we plot the 
optimized field A (£) together with the evolution of the 
populations for two particular cases, with K = 2 and 
K = 3. At first sight, it can be seen that the field shows 
oscillations wich are mounted on a step-like function. 
The latter feature is preserved from the sudden-switch 
field, wich we used as an initial guess for the optimiza¬ 
tion. Fourier transform of the driving signal reveals that 
there is only one dominant frequency / e , which together 
with the maximum amplitude A max , characterizes the 
overall shape of the field. Remarkably, this behaviour is 
common to all high-order control processes studied in 
our model, even for K > 3. In order to quantitatively 
analyze the driving field, we studied the dependance of 
f e and A max as a function of the distance eo between 
the ACs. Results are shown in Fig. [B] for K = 2, where 
the linear dependence of both quantities with e 0 is clear 
and can even be regarded as exact for the frequency, for 
which we can write f e = 

The regular behaviour shown by the numerically op¬ 
timized control field has some interesting consequences. 
First, note that as eo increases and the avoided crossings 
get further apart, the driving field will require a bigger 
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Figure 4. (color online) (a) Initial and optimized control 
fields A(t) for process P 2 , using A 1 /A 0 = 1 and eo/Ao = 10. 
Inset show the time evolution of the populations for each one 
of the diabatic states |fc) (k = 0,1, 2), given by the optimized 
field. Evolution time is set at T — Tq§ l ~ 0.91Tg 2 ^. (b) 
same as (a) for process P3, with same parameter values and 
A 2 /A 0 = 1 as well. Evolution time is set at T = Tq$ l ~ 
0.85Tg 3) . 


intensity and a larger bandwith in order to be imple¬ 
mented. In practice, at some point this requirement will 
no longer be fulfilled, and most likely the QSL time will 
tend to Tg K> for all practical purposes. This is indeed 
reasonable, since technical limitations would then imply 
that the ACs are effectively isolated, with no possible 
coupling between them. However, from a theoretical 
standpoint, this is a much different scenario than the one 
usually obtained in QOC optimizaton, where the broad 
bandwidth requirements originates from the highly 
irregular features of the optimized field. In our case 
the control function A(f) can be readily described by a 
few parameters. We associate this remarkable feature 
with the special characteristics of our model, which 
shows localized two-level interactions in a many-level 
spectrum, a scenario which is common in many different 
physical setups, as previously mentioned. 


B. Analytical approximate solution for the 
time-dependent problem 

The shape of the field also gives us interesting insight 
about the physical mechanisms involved in the observed 
enhancement of the QSL time for these control processes 
[40] , Remarkably, we found that an anallytical approx¬ 
imation for the time-dependent evolution can be drawn 
inspired from the results of the optimization process. We 
will show this solution in the following for the case K = 2, 
although the idea can extended higher order processes. 
Recall the Hamiltonian H 3 ( A) from Eq. @ , which can 
be written as the sum of its non-diagonal and diagonal 



Figure 5. (color online) Frequency f e (right axis) and max¬ 
imum oscillation amplitude A max (left axis) of the optimal 
field for process P 2 (see Fig. |4| as a function of eo- Dashed 
lines indicate lineal dependences of both quantities with eo- 


parts 


H 3 {X) = H nd + H d (X), (16) 

in such a way that IL/vd depends on the coupling param¬ 
eters Ao and Ai while the dependence on the control 
parameter is concentrated in Hjj( A). We propose the 
following expression for the driving field 


m = 


Aacos (u>t + <f>) 

eo + A^cos (u(t - t m ) + 


0 < t < t m 
t m <t<T 


(17) 

This field has the form of a step-wise constant function 
with oscillations of angular frequency ui mounted on each 
step (note that, from the previous analysis, we can infere 
that oj = eo). The held then oscillates around a fixed 
value at each step, corresponding to the localization of 
the two ACs: a t = 0 it begins at A = Ao = 0, and then 
turns to A = Ai = eo at some t = t m . The overall shape 
of A (t) then emulates the optimized held seen in Fig. 4] 
(a), with the difference that we use a constant amplitude 
Aa for the oscillating term, for convenience. 


Let us hrst consider the dynamics from t = 0 to t = 
t m . We propose that the total evolution operator for this 
evolution can be factorized as 


U 0 {t,0) = uk A) Wk B) (t), (18) 

where Uo(A)(t) = exp i f* m Hn(t')dt is diagonal in 
the diabatic basis and the superscript emphasizes the 
fact that we are working on the AC located at A = 
Ao = 0. The problem is then to hnd the unitary opera¬ 
tor U^ B \t), which satishes the Schrodinger equation in 
the interaction picture with 










































(A)t (A) 

Hnd = U 0 HmdUq being the corresponding trans¬ 
formed Hamiltonian, which takes the form 

iAi / 0 e 2iA ‘A 0 0 \ 

H' ND (t) = — A 0 0 e ie ° t A 1 , (19) 

1 \ 0 e 2l ( A *- e °‘)Ai 0 ) 

where we have defined \ = A*(t) = f* m \{t')dt' = 
^sin (ujt + (/)) — (j )o and 0o = y^sin (</>). The unitary 
evolution problem is then casted in terms of this time- 
dependent Hamiltonian. The key to consider here is that 
the exponentials that appear in the previous expression 
can be written in Fourier series using the identity 


giz sin 7 


E 


n =—oo 


( 20 ) 


where J n (z) simbolyzes the Bessel J-function of order n. 
The results we obtained from the QOC procedure indi¬ 
cates us that the frequency of the driving at each step lo 
is much larger than Ao, Ai. Then, most of the terms in 
Eq. (19) oscillate very quickly and can thus be neglected. 
This kind of rotating wave approximation is typically in¬ 
voked when analyzing high-frequency modulation of pe¬ 
riodic potentials [53] , for example in a cold atoms setup 
m ■ Formally, we approximate 


i(A;(i)-e 0 i) _ e ic /> 0 ^ f l\ e -i( 

n \ ^ / 

= e i(0o-0)j 1 A 1} 


\ n —i(nuj—eo)—in(j) 


( 21 ) 


U! ) 


where we used the argument of the previous paragraph 
to identify the term n = 1 as the resonant one and set 
u! = eo, which was expected from the numerical analysis 
of the optimal fields. Its straightforward to calculate the 
rest of the elements of Hamiltonian of Eq. (19), which 


can be approximated by a time-independent expression 


(18) where Uy B \t) = exp (— iH' ND t). Note that this 
factor introduces the couplings between the diabatic 
states which generate the time evolution of the opera¬ 
tors. The role of the driving in this process is clear. In 
the absence of the oscillatory field, i.e. A a = 0, Eq. (23) 
gives Aq = A 0 and A' : = 0, and so only states |0) and 
|1) can be connected in this evolution. This is exactly 
what we expected from the sudden-switch protocol and 
the adiabatic elimination procedure discussed in Section 
[Hi When the oscillatory field is turned on, A 0 decreases 
and Ai takes a non-zero value, thus coupling weakly 
states |1) and |2). This generates an evolution where the 
goal state of the protocol | 2 ) can draw a portion of the 
population of the other levels even when the dynamics is 
mainly dictated by the first AC. Thanks to this feature, 
the evolution towards the goal state is accelerated, thus 
providing the overall enhancement of the QSL time 
shown in the previous Section. Note that Eq. (23) 
formalizes the fact that the amplitude A a cannot be 
neglected with respect to eo- Moreover, setting Aa/co to 
a constant value is consistent with the analysis shown 
in Fig. [4] (c), where both quantites showed a linear 
correlation. 


An analogous procedure can be done for the evolution 
between t = t m and t = T, in such a way that we can 
finally write for the whole evolution as 


f U^ A \t)U^ B \t) , 0 <t<t m 
= < u[ A \t)u[ B \t)x , (24) 

l U^ A) (t m )U^ B \t m ) , t m < t < T 

where we have defined u[ A \t) = exp (—i f f T H£>(t , )dt / 'j 

and U[ B \t) = exp (— iH^ D (t — t m )) with an effective 
time-independent Hamiltonian given by 

/ 0 e-^-^A" 0 \ 

H’ nd = - e'^-^A" 0 e^°A" . (25) 

\ 0 e ~ i( Po a " 0 J 


/ 0 e-^Af, 0 \ 

H' nd = -[ e^°A ' 0 0 e^-^Ai , (22) 

\ 0 e-^-^Ai 0 J 

where we have introduced the renormalized interaction 
rates 

A ° ^ J ° (^) A ° 

Ai = J, (^j Ai (23) 

Then, the evolution of the system for 0 < t < t m is 
completely determined by the evolution operator in Eq. 


The renormalized interaction rates are now inter¬ 
changed with respect to the previous case, 

A ” = Jl (^) A ° 

A'/ = Jo ^ 7 ^^) Ar, (26) 

which is natural since in the second step the dominant 
interaction is due to the AC between states | 1 ) and | 2 ). 
In Fig. [ 6 ] we show the time evolution of the populations 
for different cases, as predicted by the analytical for¬ 
mula ( [24] ) . There, it can be seen how this expressions 
approximates very well the optimized evolution, even 
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Figure 6. (color online) Top: Optimized and analytical pro¬ 
posal of Eq. (171 for the control field as a function of time. 
Bottom: time-evolution of the populations for each of the 
diabatic states. Thin dash ed lines corresponds to the ana¬ 
lytical solution of Eq. (241. (a) Case eo = 10 Ao. (b) Case 
eo = 20 Ao. 


though the driving fields are not exactly equal. This 
behaviour allows us to assert that the high-frequency 
oscillations of the driving field at each AC, switches on 
the adjacent ACs allowing for the population of other 
energy levels and thus providing the mechanism for the 
overall speed-up of the control processes. Note that 
the solution we provide here is based on the process 
Pk= 2 ) being the next step in complexity of the two-level 
case, where this novel effects are, of course, absent. We 
believe that a similar procedure could be applied to find 
analytical expressions for higher-order processes. 


C. Initial guess and performance of the 
optimization 

We recall that the previous analytical discussion was 
motivated by the fact that the control field obtained 
from QOC had a simple shape. This, in turn, related 
with the fact that the optimized field preserved certain 
features of the initial guess Ao(i) we employed, for 
example, the step-wise structure. We will now briefly 




Iteration step m Iteration step m 


Figure 7. Top: Optimized control field as a function of time. 
Bottom: Infidelity as a function of the step number m for 
each optimization, (a) Initial guess is a linear function, (b) 
Initial guess is a sinusoidal function. 


discuss the role of the initial guess in the optimization. 
Note that the results we showed in Section hyi are 
independent of the initial guess we propose for the field. 
However, this election changes the overall performance 
of the optimization. In Fig. [7] we show the infidelity and 
optimized field obtained for two choices of Ao (f) different 
from the one shown in Fig. [4] In one of the cases, Fig. 
[7] (a), the field has a linear dependence and connects the 
positions of the ACs. The corresponding optimized field 
develops fast oscillations and overall looks very similar 
to the one in Fig. [4j Moreover, the frequency of the 
oscillations f e is the same for both cases, and the total 
number of iterations required for the convergence of 
the algorithm is also of the same order (around 4000). 
On the other hand, in Fig. 7] (b) we used as an initial 
guess a sinusoidal field, with initial and final value at 
A = A 0 = 0, and for which we deliberately change the 
parity with respect to the other cases. In that case, 
the number of iterations required by the optimization 
to converge raises by a factor of 10. Moreover, the 
optimized fields has a very irregular shape, showing 
peaks of very large amplitude (50 times bigger than the 
other cases). 

As we mentioned previously, the fact that the QSL 
time for control processes can be drawn from the opti¬ 
mization procedure enforces the power of QOC as a tool 
in this context. The results we show here also tells us 
that the performance of QOC, and its ability to give us 
information about the physical mechanisms involves in a 
control processes can be enhanced by properly providing 
the optimization with a good initial guess. In this case 
we have done so by analyzing the characteristics of the 
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system, and more precisely by studying the structure of 
the energy spectrum. 


VI. FINAL REMARKS 

In this work we studied time-optimal quantum control 
in systems with multiple avoided crossings in their energy 
spectrum. Based on previous works [38] [39], we are able 
to ellaborate protocols which generated the desired con¬ 
trol processes. We used these protocols as initial guesses 
for an optimization procedure which leads us to several 
results. Following recent work .40], we discussed how the 
quantum speed limit time for these systems as a function 
of the separation between the AC is enhanced with re¬ 
spect to the result derived from the two-level approxima¬ 
tion. We observed that when the dimension of the system 
is increased (and with it, the number of ACs involved in 
the dynamical process) the speed-up becomes more pro¬ 
nounced. Having obtained the control protocols at the 
QSL, we numerically analyzed the shape of the fields de¬ 


rived from the optimization. We found that they showed 
a very regular behaviour, characterized by the presence 
of single-frequency oscillations mounted on a step-wise 
function. Based on these results, we were able to con¬ 
struct a model for the control problem which we solved 
analytically. This model also allows us to qualitatively 
explain that the main feature behind the optimization 
was the collective dynamics of multiple avoided-crossings. 
Finally, we studied how the outcome and performance of 
the optimization were modified when varying the initial 
guess for the control protocol. We found that using dif¬ 
ferent initial guesses can lead to very different shapes of 
the optimized control field. This behaviour indicates that 
a preliminary analysis of the system spectrum, as done 
here, can act as pre-optimization method as it lead us to 
a good choice of the initial guess. 
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